Import libraries

In [1]:
from glob import glob  
from natsort import natsorted, ns

import geopy
import folium
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter

import rasterio  
from rasterio.plot import show
from pyproj import Transformer
from rasterio.windows import Window
import rioxarray 

import imageio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

%matplotlib inline

Search all tifs and sort files

In [2]:
def tif_file(path):
    
    tif_file =[]
    
    # using glob library to get all the file with .tif
    files = glob(path,recursive = True) 
    for file in files: 
        tif_file.append(file)
        
    # sort files with number in the file (ascending order) using natsorted library
    tif_file = natsorted(tif_file, alg=ns.IGNORECASE)     
    return tif_file
In [3]:
# pwd '/Users/adamtky/Desktop/BeCode/3. Project/Project 3 - 3D Houses'
In [4]:
path = './DSM/**/*.tif'
DSM_tif = tif_file(path)
DSM_tif 
Out[4]:
['./DSM/DSM PROV_BRABANT_WALLON/Brabant_Wallon_dsm_1.tif',
 './DSM/DSM_PROV_HAINAUT/Hainaut_dsm_2.tif']
In [5]:
path = './DTM/**/*.tif'
DTM_tif = tif_file(path)
DTM_tif 
Out[5]:
['./DTM/DTM PROV_BRABANT_WALLON/Brabant_Wallon_dtm_1.tif',
 './DTM/DTM_PROV_HAINAUT/Hainaut_dtm_2.tif']

Geocoding - single address

In [6]:
address = 'Abbaye de Bonne-Espérance'
In [7]:
# enter an address in Belgium
# address = input("Enter an address in Walloon Brabant or Hainaut , Belgium:")

def lat_long(address):
    
    # to get the longtitude and latitude of the address entered
    locator = Nominatim(user_agent="myGeocoder")
    location = locator.geocode(address)
    house_lat_long = [location.latitude, location.longitude]
    return house_lat_long

def house_map(func):
    
    # to plot address on a map
    house_map = folium.Map(location=func,zoom_start=150)
    folium.Marker(location=func,popup=lat_long(address)).add_to(house_map)
    house_map.save("house_map.html")
    return house_map 

house_map(lat_long(address))

# Example address: 
# Abbaye de Bonne-Espérance
 
Out[7]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Convert longtitude & latitude

In [8]:
# longtitude and latitude of the address entered
lat,lon = lat_long(address)
lat,lon
Out[8]:
(50.38571435, 4.144005338377226)
In [9]:
def transform(lon, lat):
    
    # transform to Belgium 'EPSG:31370' coordinate
    transformer = Transformer.from_crs("EPSG:4326", crs_to = 'EPSG:31370' ,always_xy=True)
    xx, yy = transformer.transform(lon, lat)
    print(f'longtitude:{xx} latitude:{yy}' )
    
    return xx,yy
In [10]:
# coordinate EPSG:31370
xx,yy = transform(lon, lat)
longtitude:134017.30872732686 latitude:119445.92955795769

Locate tif for given longtitude & latitude

In [11]:
# create all bounding box from tifs

def bounds(tifs):
    
    bounds = []
    
    for i in tifs:
        src = rasterio.open(i)
        bounds.append(src.bounds)
        
    return  bounds 

src_bounds = bounds(DSM_tif)
In [12]:
src_bounds[:3]
Out[12]:
[BoundingBox(left=130000.00000000211, bottom=133999.99999999837, right=196000.0000000021, top=167999.99999999837),
 BoundingBox(left=42000.0, bottom=69999.99999999927, right=168000.0, top=167999.99999999927)]
In [13]:
scr_path = []

if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):
    
    scr_path.append(DSM_tif[0])
    
    print(f'The house is in this tif : {DSM_tif[0]}') 
    
elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):
    
    print(f'The house is in this tif : {DSM_tif[1]}') 
                                                                                                       
else:
    print('The house is not found in the tif')                                                                                                    
The house is in this tif : ./DSM/DSM_PROV_HAINAUT/Hainaut_dsm_2.tif
In [14]:
# check which tif contains the location in DSM

def check_DSM(xx,yy):
    
    scr_path = []

    if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):

        scr_path.append(DSM_tif[0])

        print(f'The house is in this tif : {DSM_tif[0]}') 

    elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):

        print(f'The house is in this tif : {DSM_tif[1]}') 
        
        scr_path.append(DSM_tif[1])

    else:
        print('The house is not found in the tif')    
            
    return scr_path

dsm_path = check_DSM(xx,yy)
The house is in this tif : ./DSM/DSM_PROV_HAINAUT/Hainaut_dsm_2.tif
In [15]:
# check which tif contains the location in DTM


def check_DTM(xx,yy):

    scr_path = []

    if (xx >= src_bounds[0][0] and xx <= src_bounds[0][2]) & (yy >= src_bounds[0][1] and yy <= src_bounds[0][3]):

        scr_path.append(DTM_tif[0])

        print(f'The house is in this tif : {DTM_tif[0]}') 

    elif (xx >= src_bounds[1][0] and xx <= src_bounds[1][2]) & (yy >= src_bounds[1][1] and yy <= src_bounds[1][3]):

        print(f'The house is in this tif : {DTM_tif[1]}') 
        
        scr_path.append(DTM_tif[1])

    else:
        print('The house is not found in the tif')    
            
    return scr_path

dtm_path = check_DTM(xx,yy)
The house is in this tif : ./DTM/DTM_PROV_HAINAUT/Hainaut_dtm_2.tif

Clip the location of the house from tif

In [16]:
dsm_path
Out[16]:
['./DSM/DSM_PROV_HAINAUT/Hainaut_dsm_2.tif']
In [17]:
dtm_path
Out[17]:
['./DTM/DTM_PROV_HAINAUT/Hainaut_dtm_2.tif']
In [18]:
def clip_dsm(path,window_size:int):
 
    xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
    
    # set window size
    n = window_size

    # create coordinates and geometries
    coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]

    coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]

    geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]

    # clip the image as per the geometries size
    clipped = xds.rio.clip(geometries)

    clip_dsm = clipped.rio.to_raster("clipped_dsm.tif",tiled=True, dtype="int32")

    return clipped.plot()
In [19]:
def clip_dtm(path,window_size:int):
 
    xds = rioxarray.open_rasterio(path,masked=True,chunks=True)
    
    # set window size
    n = window_size

    # create coordinates and geometries
    coor1,coor2 = [(xx-n),(yy+n)],[(xx+n),(yy+n)]

    coor3,coor4 = [(xx+n),(yy-n)] ,[(xx-n),(yy-n)]

    geometries = [ {'type': 'Polygon', 'coordinates': [[coor1,coor2, coor3,coor4,coor1 ]]}]

    # clip the image as per the geometries size
    clipped = xds.rio.clip(geometries)

    clip_dtm = clipped.rio.to_raster("clipped_dtm.tif",tiled=True, dtype="int32")

    return clipped.plot()
In [20]:
dsm_path[0] 
Out[20]:
'./DSM/DSM_PROV_HAINAUT/Hainaut_dsm_2.tif'
In [21]:
dtm_path[0]
Out[21]:
'./DTM/DTM_PROV_HAINAUT/Hainaut_dtm_2.tif'
In [22]:
clip_dsm(dsm_path[0],50)
Out[22]:
<matplotlib.collections.QuadMesh at 0x7f7eebc0ca50>
In [23]:
clip_dtm(dtm_path[0],50)
Out[23]:
<matplotlib.collections.QuadMesh at 0x7f879032f150>
In [24]:
def chm_tif():
    
    # open the digital terrain model 
    
    with rasterio.open('clipped_dtm.tif') as src:
        lidar_dtm_im = src.read(1, masked=True)
        dtm_meta = src.profile
        
    # open the digital surface model 
    
    with rasterio.open('clipped_dsm.tif') as src:
        lidar_dsm_im = src.read(1, masked=True)
        dsm_meta = src.profile
    
    # calculate canopy height model
    
    lidar_chm = lidar_dsm_im - lidar_dtm_im
    
    with rasterio.open('clipped_chm.tif', 'w', **dsm_meta) as ff:
        ff.write(lidar_chm,1)
        
    chm_tif = 'clipped_chm.tif'
    
    return chm_tif

Create 3D

In [25]:
def House_3D(tif,size:int,azim=215):

    chm = imageio.imread(tif)

    nx,ny = chm.shape

    x = np.linspace(0, size*2, nx)
    y = np.linspace(0, size*2, ny)

    yv,xv = np.meshgrid(x, y)

    fig = plt.figure(figsize=(15,15))
    ax = fig.add_subplot(111, projection='3d')

    chm3d=ax.plot_surface(xv,yv,chm,cmap='plasma',linewidth=0)
    ax.set_title('CHM(Canopy Height Model)')
    ax.set_xlabel('Distance (m)')
    ax.set_ylabel('Distance (m)')
    ax.set_zlabel('Elevation (m)')
    ax.view_init(azim=azim)
    
    fig.colorbar(chm3d, shrink=0.3, aspect=10);
    fig.savefig('house.png', dpi=200) 
    
    return  plt.show()
In [26]:
House_3D(chm_tif(),50,300)
In [27]:
# compare google map image of the house

from IPython.display import Image
Image("Basilique Notre-Dame de Bonne-Espérance.png")
Out[27]:
In [28]:
# to view actual house on google map

import webbrowser
url = 'https://www.google.com.my/maps/place/'+str(lat)+','+str(lon)
webbrowser.open(url)
Out[28]:
True
In [29]:
from osgeo import gdal
from mayavi import mlab

ds = gdal.Open(chm_tif())
data = ds.ReadAsArray()

data = data.astype(np.float32)

mlab.figure(size=(640, 800), bgcolor=(0.16, 0.28, 0.46))

surf = mlab.surf(data,warp_scale='auto') 

#mlab.show()
In [30]:
Image("Basilique Notre-Dame de Bonne-Espérance_mayavi.png")
Out[30]: